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Abstract. The bulk macroscopic response of a system of particles or inclusions with 
field-induced forces is studied. The susceptibilities and transport coefficients in such 
a system are expressed as averages of a multiple scattering expansion. A special 
diagrammatic method is developed to analyze the structure of the expansion. The 
concept of irreducibility is discussed in detail and shown to be crucial in obtaining 
macroscopic equations characterizing the system response with coefficients depending 
solely on local properties of the medium. Due to the representation of particles by lines 
in diagrams, irreducibility is given a particularly simple topological interpretation in 
the diagrammatic language. The method is illustrated by a discussion of response 
problems in colloidal suspensions in presence of hydrodynamic interactions. 



PACS numbers: 66.00.00,47.57.-8 



Diagrammatic approach to response problems in composite systems 



2 



1. Introduction 

Calculation of the effective properties of a composite system is an active field of research 
(see e.g. lU [21 E] and references therein) important not only for the physical insight it 
provides but also for many potential practical applications. The composite materials 
considered here are systems of particles or inclusions embedded in a homogeneous 
medium and subject to an external field. A classical example of such a system is the 
Kirkwood-Yvon dielectric [H, [5] - a set of polarizable, spherical inclusions embedded in 
a uniform and isotropic medium. Its relative simplicity makes it a convenient starting 
point to illustrate the methods presented here. Next, we focus on a more complicated 
composite system - a colloidal suspension, in which the motion of suspended particles 
in the liquid is caused either by the gravitational force or by an imposed external flow. 

To define the effective macroscopic properties, one must start from local equations 
that govern the system response to external disturbances. The construction of such 
equations is not trivial if there are long-range interactions present, since they often 
lead to the divergent integrals in the expressions for the transport coefficients. Those 
divergences are usually removed with use of rather subtle "regularization" techniques 
(e.g. [HI [71 El E] ) , which involve nontrivial manipulation of multiple scattering expansion 
with the careful resummation of the various kinds of terms. The calculation may be 
facilitated by the development of diagrammatic methods which not only allow the local 
response equations to be obtained quickly and reliably, but also provide us with a clear 
interpretation of the different steps in the regularization procedure, which are sometimes 
obscured in the standard approach. 

A key factor for a successful diagrammatic method is the requirement that the 
structure of the terms of the scattering expansion should be reflected in topological 
properties of respective diagrams. In particular, since a given particle may take part 
in more than one scattering event, it is convenient to represent particles not by points 
but by lines in analogy to the diagrammatic techniques developed by the Brussels group 
[To] in nonequilibrium statistical physics. In particular, due to the representation of the 
particles by lines in our diagrammatic approach, a natural ordering of the successive 
scattering events in the multiple scattering expansion is reflected in the ordering of 
the scattering events along the particle line. Additionally, the notion of irreducibility, 
central to the regularization procedure, is now given an elegant interpretation in terms 
of the topology of the diagrams. This constitutes a fundamental difference between 
our approach and another diagrammatic technique found in the literature, due to 
Barrera [TTl [T2] . In Barrera approach the particles are represented by points, which 
complicates the analysis, since the diagrams then become multiply- connected, i.e. 
there is usually more than one edge linking the nodes. Hence the edges must be 
numbered in order to obtain a unique identification for a particular diagram. This 
makes it harder to analyze various types of diagrams and to link the structure of the 
multiple scattering expansion to the their topological properties. Additionally, there is 
no obvious generalization of that technique to the time-dependent case, in contrast to 
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the particle line approach. 

The regularization procedure with use of the diagrammatic technique allows one to 
obtain well-defined theoretical expressions for transport coefficients, free of the integrals 
diverging with the size of the system, even in the presence of long-range forces. In 
particular, as it will be shown in a subsequent paper, the diagrammatic expansion 
allowed us to construct a hierarchy of equations for the correlation functions in a 
settling suspension, which in turn allows to solve the long-standing problem of velocity 
fiuctuations in non-Brownian suspension [13] . Namely, it was argued theoretically more 
than 20 years ago by Cafiisch and Luke [H] that the velocity fiuctuations should diverge 
linearly with the macroscopic size of the system. However, this prediction has not 
been confirmed in the experiments [151 [IHl [H] • Instead, in most of the experiments, 
the saturation of the velocity fiuctuations was observed. A careful analysis of the 
correlation structure of the settling suspension, in which the diagramatic analysis plays 
a fundamental role, allowed us to show that the velocity fiuctuations do not diverge with 
increasing container dimensions. Another problem of a similar nature is the calculation 
of the mean velocity of a settling non-Brownian suspension. Batchelor [T8l[T9] calculated 
this quantity for the polydisperse suspension. It turns out, however, that his theory gives 
ambiguous results for the monodisperse case (the result depends on the way the limit 
is taken) [20]. Also in this diagrammatic analysis allows us to derive a well- 

defined and unambiguous result for the sedimentation velocity in both polydisperse and 
monodisperse case [13] . 

The diagrammatic expansion constitutes also a good starting point for the 
construction of various approximation methods for calculating the effective properties 
of the medium. In general, transport coefficients have different values in the short-time 
regime i.e., for times in which particles have hardly moved and for long times when the 
relaxation of the distribution of particle positions becomes important. This relaxation 
gives rise to the memory effects, which can also be incorporated into the presented 
diagrammatic approach. Additionally, we discuss the relation of our approach to another 
method of obtaining the transport coefficients, based on the Fourier space formulation 
of response equations and subsequent calculation of the small wavenumber, k — 0, limit 
of the response kernels. 

2. Multiple scattering expansion 

A composite medium is often modeled by a disordered system of particles or inclusions 
embedded in a homogeneous matrix. In many cases, if such a system is inserted into 
the field \I/o(^), the particles themselves become sources of the field (as it is the case for 
polarizable dipole systems). The contribution of the induced sources to the total field 
in the sample, \l/(r) is then given by 




(2.1) 
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where the function s{r) describes the intensity of the sources and G{r, r') is the Green's 
function. The response of the particle i to the field "^ext, external to the particle, is 
characterized by the operator M 

s,(r) = J dr'M{i; r, r')^e.tir'), i = 1, . . . , N (2.2) 

with 

M{t; r, r') = ^,(r)M(2; r, r')e,{r') (2.3) 

where 6i is a characteristic function of zth particle. The above equation reflects the fact 
that the induced sources Si depend only on the values of the field "ifextir') within the 
particle i and vanish outside the particle. 

These ideas may be illustrated with an example of the Kirkwood-Yvon dielectric 
[H E] - a system of N identical polarizable point dipoles. In this case the field ^ 
corresponds to the electric field in the dielectric whereas the sources Sj(r) are expressed 
in terms of the dipole moments, pi, as 

s,(r) = p,5(r - r,) (2.4) 

The Green's function is then given by dipole-dipole interaction tensor 

G{r,r') = g{r-r'), 

, „„1 1 3ff (2-5) 
6? r = VV- = -- + ^. 

Finally, the single particle scattering operator is simply 

M{i-r,r') = 5{r - Ri)a5{r' - Ri) (2.6) 

where a is the molecular polarizability and Ri - position of i th dipole. 

Let us find the response of a composite system to the external field '^Q{r). The 
total field in the sample is then given by 

m = m^ + Y,Gsi (2.7) 

i 

whereas 



^o + GY^sA 



Si = M{t){^o + G}_^s,] {2.t 

In the above, the shorthand notation is used, in which the integrations and the 
coordinates (r, r') are suppressed, i.e. 

{AB){r) = j A{r,r')B{r')dr' (2.9) 

Additionally, the space arguments of the operators (r,r' etc.) are dropped. Note that 
the term i = j is omitted in the summation (12.81) since the response relation (12. 2p relates 
the sources Sj to the field external with respect to the particle i. 
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The relation (12.81) is iterated to obtain successive terms of a multiple-scattering 
expansion 

Si = (M(2) + J2 M{t)GM{j) + ... )^o. (2.10) 

Using the above formalism, one can find the response kernel T defined by the relation 

s=T^o, (2.11) 

where s is the total source intensity 

s = J2si- (2.12) 

i 

Using fl2.10p one may represent T in the form of a scattering expansion 

T = Y, M{{) + E E ^(OGM(j) + . . . (2.13) 

i i j^i 

3. Averaging the scattering expansion over configurations 

Since we are interested in the average behaviour of the system on a macroscopic level, 
its response should be averaged over an ensemble of particle configurations. Averaging 
of (12J31) leads to 

< s >=< T > ^0 (3.1) 
where the brackets stand for a configurational average 

A>= j A{X,r,r')P{X)dX, (3.2) 



< 



and P{X) is the configurational probability distribution function, with X = 
{Ri, R2, • • . , Rn}- 

In the dielectric example considered above this corresponds to the relation between 
the external electric field, Eq, and the polarization, P, 

P =< s >=< T > Eq. (3.3) 

However, the above relation is not local, since polarization in the sample depends not 
only on Eq but also on the shape of the sample, boundary conditions etc. Conversely, a 
local relation characterizing dielectric response is 

P = eox < E > (3.4) 

where < E > is the macroscopic electric field. The electric susceptibility x does not 
depend on the shape or size of the sample but only on the local properties of the material. 
In particular, the dielectric constant of a medium is expressed as 



e = 1 + X 



(3.5) 
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Going back to the general case, we see that the operator < T > may not be a physically 
meaningful measure of system's response. Instead, one should study the response 
operator X defined by the relation 

<s>=X<^> (3.6) 

linking the sources < s > to the total field inside the sample, < \E' >. The procedure of 
obtaining X from < T > (so-called "reduction" or "regularization" of a response kernel) 
is presented below in a systematic way. 
First, we represent the operator T as 

^ = E^(^) + ^E^(^'^') + ^ E n^,j,k) + ..., (3.7) 

where T(zi, ...,««) comprises all these terms in the scattering sequence in which all the 
particles {zi,Z2---^s} are included. 
Hence we get for < T > 

^ ATI /■ 

< T >= J] y T(l, 2..., s)P(l, . . . , iV)dl...diV, (3.8) 

where we write i instead of Ri to simplify notation. The above expression may also be 
written as 

TV 

<T>=J2- Til, 2..., s)n(l, 2, . . . , s)dl...ds, (3.9) 

8 = 1 

where n(l,2...,s) is the s-particle partial distribution function 

A^! f 

n(l, 2, s) = J P(l, .., iV)d(s + l)...dN. (3.10) 

Note that the s-particle partial distribution function can be written as 

n(ri,r2,...r,) =< J^' 6{r, ~ R,,)6{r2 - R^,)...6{r, - R^J >, (3.11) 

11,12, ■■■As 

which in a shorthand notation will be also denoted as < 1 2 ... s >. The sum in 
the above expression is supplied with the condition that all different 
each from the other. The above definition fl3.1ip of partial distribution function holds 
also for a system with a variable number of particles (if the grand canonical ensemble 
is used). In this case the sum in (13.91) should be extended up to infinity: Xli^i W^ ■ 

Next we assume that the correlations between the two groups of particles vanish 
as the distance between them goes to infinity. This means that the partial distribution 
function should have the group property, i.e. 



n{l, 2...,r, r + 1, .., s) — > n{l, 2...,r)n{r + 1, s). 



(3.12) 
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as the distance between the particles {1, 2...r} and {r + 1, s} goes to infinity. 

This property of the partial distribution function allows us to decompose n(l, 2...s) 

as [in! 

n(l) = /i(l), 

n(l,2) = n(l)n(2) + /i(l,2), 

n(l, 2, 3) = n(l)n(2)n(3) + n{l)h{2, 3) + n(2)/i(l, 3) + n(3)/i(l, 2) + 2, 3), 

(3.13) 

where the s-particle correlation function h{l,2...s) which vanishes whenever any subset 
of particles C {1, 2...s} is dragged away from the rest. 

With the above decompositions one can write the average < A > in ( \3.2\\ as the 
sum of terms of the general form 

Ts{A, c) = j Z2, • • • , is)c{ii, 12, ■ ■ ■ , is)dii...dis, (3.14) 

where c{ii,i2, ■ ■ ■ ,is) is a product of a number of correlation functions involving 
particles {ii, . . . ,is} whereas A(ii, ^2, • • • , is) is one of the scattering sequences making 
up T{ii,i2, . . .,is)- 

For the later use we introduce after Michels p2] the "uncorrelating operator" 

Punc=><, (3.15) 

which has the property of statistically uncorrelating the variables at its left from those 
at its right, i.e. 

< APuncB >=< A><B> . (3.16) 
The orthogonal complement of Punc is 

Qunc = !-><. (3.17) 

So, using the notation of Eq. (13. lip , we get for example 

< 1 Qunc 2 >=< 1 2 > - < 1 >< 2 >= n(l, 2) - n{l)n{2) = h{l, 2). (3.18) 

The decomposition (13.131) together with the cluster expansion (13.71) leads to the 
representation of response kernels as sums of many-body terms from the scattering 
sequence multiplied by respective correlation functions. To deal effectively with such a 
complicated structure a special diagrammatic technique is employed. 

4. Diagrammatic representation 

We introduce the diagrammatic representation of the scattering (S) and correlation (C) 
structure of the kernels. Such SC diagrams consist of the following elements 



Diagrammatic approach to response problems in composite systems 



8 



(i) the horizontal hne represents a given particle (also called particle line) 

(ii) the symbol O stands for the operator M{i) 

(iii) the vertical hne | stands for the G - bond 

(iv) double vertical line represents the correlation function h (called h-bond) 

The exact interpretation of an h-bond depends on the geometric structure of a 
diagram. For example 
lo 

2) 

3I 

stands for h{l, 2, 3), whereas 
1( 



corresponds to /i(l, 2)/i(2, 3). 

Moreover, if the first symbol on the particle line (looking form the left side) is filled, 
then the position of this particle is integrated over. Hence, for example the diagram 



10 



(D 1) 



represents the expression 



/ 



dld3d4/i(13)M(l)G(12)M(2)G(23)M(3)G(31)M(l)G(14)M(4). 



(4.1) 



Note that the diagrams should be read from left to right. The particles lines 2 and 
4 in the above diagram are left out since there's only a single operator involving each of 
these particles. 



4.1. Irreducibility 

A key notion in the analysis of internal structure of scattering sequence terms is the 
concept of irreducibility of a diagram. Namely, the G bond in the diagram is called a 
connection line if the removal of this G-bond causes the diagram to become disconnected. 
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Diagrams with one or more connection lines are called reducible, whereas diagrams 
without any connection lines - irreducible. 
For example the diagram 



b 



(D 2) 



is reducible and can be cut into two pieces by breaking the G - bond between particles 
2 and 3. The connection line which is most to the left will be called articulation line. 
Note that the sub-diagram on the left of the articulation line is irreducible. 

In the analogous way we can define the reducibility for the scattering structure of 
the diagrams (S-reducibility). First of all nodal line is defined as a G-bond which 
would be a connection line if all the h — bonds in a diagram are removed. Diagrams 
with one or more nodal lines arc called S-reducible. 

Hence in the following diagram 



(D 3) 



the G bond between particles 2 and 3 is the nodal line but not the connection line and 
the diagram is S-reducible, although it is irreducible with respect to its full SC-structure 
(which includes both correlation and scattering part). 



4-2. The nodal structure 

The nodal lines decompose the particles in a given diagram on the set of nodal blocks 
Cf. Ci denotes the set of particles on the left of the first nodal line, C2 - the particles 
between the first and the second nodal line and so on. Note that the definition of the 
nodal fine assures that Cj fl Cj = if only i j. 
For example the diagram 
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b 



d 



[D 4) 



has the nodal structure of the form 





2,3 




1 



or simply 1|23. 

The structure in the above figure is called the nodal structure graph (NSG). The 
vertices of such a graph are nodal blocks, whereas the bonds in this graph are created 
by nodal lines. 



4-3. The block distribution function 

Consider all the irreducible diagrams which have the same scattering structure and differ 
only in correlation structure. The task of summing all of these diagrams thus boils down 
to finding the sum of all their correlation functions. 

To start with, the condition of irreducibility requires that if there is a nodal line in 
the diagram then particles on the left of it cannot be totally uncorrelated from particles 
on its right. This means that the correlation function that we are looking for is given 
by 



b{Ci\C2\--.\Ck) —< Ci(l — Punc)C2{l — Pu„c)...(l — PuncjCk > 



(4.2) 



Here Ci\C2\..-\Ck describes the nodal structure of the diagram, whereas the operator 
Punc is the " uncorrelating operator" introduced in f l3.15p . The function b{Ci\C2\---\Ck) 
defined in (14.21) is called the block distribution function [7]. Note that if there are no 
nodal lines in the scattering structure of a given s-particle diagram, than b would be 
just the full s-particle partial distribution function n(l, 2, s). 



Diagrammatic approach to response problems in composite systems 11 

To get a better grip on b{Ci\...\Ck), let us evaluate it for a few simple scattering 
sequences. For the sequence presented in Diagram (D 4) the block distribution reads: 

6(1|23) =< 1(1 - P„nc)23 >=< 123 > - < 1 >< 23 >= n(l, 2, 3) - n(l)n(23). (4.3) 

We see that 6(1|23) goes to zero as the particle 1 is dragged away from the particles 2 
and 3, as in this case 

n(l,2,3) ^ n(l)n(23). (4.4) 
Let us consider now the scattering sequence of the form 





d 



where 



«i,«2 ■ --ik 



stands for any irreducible scattering sequence that involves the particles ii, 22 • • • ik- 

The above scattering sequence has the nodal structure (1|23|45). Therefore its 
block distribution function reads 

6(1|23|45) =< 1(1 - P„nc)23(l - P„„e)45 >= (4.5) 
=< 12345 > - < 1 >< 2345 > - < 123 >< 45 > + < 1 >< 23 >< 45 >= 
= n{l, 2, 3, 4, 5) - n(l)n(2, 3, 4, 5) - n{l, 2, 3)n(4, 5) + n(l)n(2, 3)n(4, 5), 

which, as can be easily proved, vanishes whenever the particle {1} is separated from the 
rest or the group {1, 2, 3} is dragged away from {4, 5}. 
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5. Reduction of the diagrams 

In Section II we have obtained the representation of the response kernel < T > as the 
sum of terms of the form fl3.14p . Each such term may be represented as a diagram, 
according to the rules formulated above. Next, those diagrams may be divided into two 
groups: reducible and irreducible ones. Thus 

< T >=< T + <D> 

where < T^^^ > is the sum of all irreducible diagrams of < T > whereas < -D > - is the 
sum of the reducible ones. However, each reducible diagram may be written in form of 
a product: 

D = IGR (5.1) 

where D stands for the diagram under consideration, / is its part to the left of the 
articulation line and R is the part to the right of the articulation line. As follows from 
the definition of irreducibility, the diagram corresponding to I must be irreducible, since 
it does not contain an articulation line itself. For example, the diagram (D 2) is divided 
in a following way 



1 



d-b 



Here I is given by the diagram 



-o 



whereas R is given by 
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The scattering structure of both / and R diagrams is exactly the same as the 
scattering structure of the original < T > diagrams. However, due to the irreducibility 
restriction, the correlation structure of / diagrams is different: the correlation function 
which multiplies a sum of all I diagrams with the given scattering structure is given 
by the block correlation function b{Ci\...\Ck) defined in fl4.2p . Thus the sum of all R 
diagrams is just < T >, whereas the sum of all / diagrams is < T >*''^. These arguments 
lead to 

< T >=< T >"■'■ + < T >'''^ G<T> . (5.2) 

which becomes exact in a thermodynamic limit [23]. Applying both sides of the above 
equation to < \Efo > and using (12.111) one gets 

< s >=< T > ^0 =< T < T G<T>'^o (5.3) 

This equation can be combined with the average of (12. 7p 

< ^ >= ^0 + G < s >= ^0 + G < T > ^0 (5.4) 

leading to 

< s >=< T >^"^ ^0+ < T >^"'^ (< ^ > -^o) =< T >^"'^< ^ > (5.5) 

which links the sources < s > to the local field inside the sample, < \E' >. Thus the X 
operator in Eq. (13.61) may be identified with < T >*'"^. 

In the following, we consider a more general form of a response kernel, namely 

A = J2 Moii) + E E M4t)GM^{j) + E E E M^i^GM{j)GM^{k) + ... 

i i j^i i j^i k^j 

(5.6) 

which differs from (12.131) in that it contains the opening operator M^{i), the closing 
operator M^{i) and the single-particle operator Mq, which in general are different from 
M{i). 

The reduction procedure for < A > is similar to the one presented above. However, 
due to the presence of M"< (i) and M"> (i) in the scattering sequence of A, the reduction 
formula is slightly more complex than (15.21) : 

< A >=< A + <A< G<A> > (5.7) 
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where the operators and have scattering sequences 

A^ = J2 ^CO+E E M^{z)GM{j)+J2 E E M^{t)GM{j)GM{k)+. . . (5.8) 

and 

/l> = J] M{z)+J2 E ^(0<^M>(j)+ J] J] J2 M{^GM{j)GM^{k)+. . . (5.9) 
respectively. 

As an example of a response problem described by a general structure (15.61) we 
consider a colloidal suspension - a system of solid particles immersed in a fluid. 



6. Transport phenomena in colloidal suspensions 



The system under consideration consists of N identical spherical particles of radius a 
immersed in an incompressible fluid of shear viscosity rj. The particle Reynolds number 
is assumed to be small so that the inertial effects are negligible and the fluid can be 
described by Stokes equations. The sources Sj are then the force density exerted on the 
fluid by the particles whereas the role of the field \l/ is played by the fluid velocity field, 
v{r). 

As it was shown by Mazur and Bedeaux [2l] if the particles are impenetrable to 
the flow and the stick boundary conditions at their surfaces are assumed, then validity 
of Stokes equations may be formally extended inside the particles: 



r^V^v -Vp + /o(r) + /(r) = 0, 
V ■v = 0, 

v{r) = u,{r) =Ui + n,x {r- R,) 
p{r) = 



for 
for 



\r 
\r 



Ri 

R; 



< a, 

< a. 



(6.1) 
(6.2) 
(6.3) 
(6.4) 



Here /o(''") is an external force density applied to the fluid, such as gravity. Next, f{r) 



is an induced force density localized on the particle surfaces |24l EH] and Ui and fli are 
translational and rotational velocities of the particles. 

The solution of hydrodynamic equations (16. ip . (16.21) can be written as 



v{r)=vo{r) + j G{r,r')-f{r')dr', 



(6.5) 



where VQ{r) is the flow in absence of the particles and G{r, r') is the Green tensor. For 
an unbounded fluid G(r, r') is given by the Oseen tensor Gq 



G{ry) = G,{r-r') 
_ . . 1 1 + ff 
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The response of a single particle to the fluid field is described by the one-particle 
friction kernel Zo{i) 

f,ir) = J Z,(z; r, r')iu,ir') - i;,(r'))dr' (6.7) 

where Va{r) is the flow field external to particle i. The above equation is a counterpart 
of the relation (12. 2p . with the operator M corresponding to —Zo- The explicit form 
of Zo{i) for variety of boundary conditions may be found e.g. in [26]. Next, we may 
proceed in several ways. 

In a friction problem, one looks for the forces induced on the particles for the 
given flow field. This leads to the relation 



fir) = J2 W = / O • ivir') - Voir'))dr', 



(6. 



where the friction kernel 2{r, r') can be represented in form of the scattering expansion 

Z = Y^ Zoi^ -EE Zo{i)GZoi3) + ... (6.9) 

i i jy^i 

The above relations are analogous to (12. lip and (I2.13P respectively. When deriving 
Eq. (16. Sp . we used the fact that the operators Zo{i;r,r') are localized inside the 
corresponding particles, together with the condition (16. 3p . Additionally, the notation 
may be simplified further by introducing the operators Zo and Q: 

Zo,, = Z^m, g,, = G(zj)(l - 6,,) (6.10) 

which are the NxN operator matrices in the particle indices. In the above, G{ij) denotes 
the operator G placed between Zo{i) and Zo{j) in the scattering expansion (16. 9p . Here 
and below we use the script letters {Zo, Q, ^ . . .) for objects acting in the particle 
index space. With the above notation (16. 9p takes form 

z = Zo{i + gZor' (6.11) 

The above allows us to find the friction matrix <^ which is defined by the relation between 
the forces and torques acting on the particles and their velocities (in the absence of 
external flow) 

:f = cu, (6.12) 

Here J- = [7- , T) is the 6 N- dimensional vector of forces and torques acting on each of 
particles: (JF, T) = (Fi, Fa, F^v, Ti, T^v) whereas W = {U,n) is the vector 
of translational and rotational velocities of the particles tl = (t/i, Un, ^i, ^n)- 
The friction matrix, ^, may be similarly decomposed as 



c 



Diagrammatic approach to response problems in composite systems 



16 



The matrices C^'' {p,(l = t or r) are the 3Nx3N Cartesian tensors, and the superscripts 
* and ^ correspond to the translational and the rotational components, respectively. 

Subsequent analysis is facilitated by introduction of multipole expansion. Namely, 
one represents the force densities and the velocity field around ith particle as the (infinite 
dimensional) vectors of successive multipoles: 



and 



udr - Voir 



V ••• / 



/ t/, - voiRi) \ 
a - uj{R,) 



(6.13) 



V 



(6.14) 



In the above, force multipoles are obtained by the following integrations of f{r) 



fir)9,{r)dr 

[r - Ri) X f{r)ei{r)dr, 



(6.15) 



{r-R,)f{r) e,{r) 



where 



eAr) 



(6.16) 



is the characteristic function of the particle i and the overbar stands for the symmetric 
and traceless part of the tensor. 

On the other hand, velocity multipoles are obtained by the following 
differentiations: ^ 

u;{Ri) = -(V X vo)r=B, (6.17) 



1 



In the multipole notation, the operators Go and Zo become matrices. The friction 
matrix, defined in fl6.12p . relates the two lowest velocity multipoles to the two lowest 
force multipoles. Therefore it can be obtained from the multipole matrix Z by the 
following projection 

C = VZV. (6.18) 

where 7-* = {'P'^,T'^) are the projection operators extracting the two lowest moments 
from the velocity (or force) distribution, i.e. 



:F = vf, 



(6.19) 
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and 

U = Vv (6.20) 

Subsequently, we will also use the operator 7-*^ which gives the third multipole of the 
force field, i.e. 

S. = VU^. (6.21) 

and similarly for the velocity field 

g. = Vfvo. (6.22) 

Let us now find forces acting on particles in the presence of the ambient flow Vq. 
From Eq. (16.81) one gets in this case 

= CU-VZvo. (6.23) 

The above formalism can also be used to solve the mobility problem: finding velocities 
of the particles lA for given forces ^ and flow t^o. In this case, the relation (I6.23P gives 

U = C^7' + C^VZvq = ^i^ + Cvo, (6.24) 

which defines the mobility matrix 

fi = (6.25) 

together with the convection kernel C 

C = nVZ. (6.26) 

The mobility matrix, /j,, allows us to find translational and rotational velocities of 
particles in terms of forces and torques acting on them in the absence of an external 
flow 



ft tr 
fl fj, 



Finally, let us consider a problem of finding the force density / for given forces 
^ 7^ and ambient flow Vq. In this case, from (16.231) and (16. 8p we obtain 

f = CT-ZvQ. (6.28) 

where C is the transpose of C operator 

C = ZVfi, (6.29) 

while the convective friction kernel Z [271 given by 

Z = Z- ZVuVZ. (6.30) 
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The operator G2 produces the velocity fields, which are force-free and torque- free. 
The scattering expansion for the convective friction kernel Z is found to be 

oo 

Z = Zo{l + gZoV = Yl ^o{-QZa)\ (6.31) 

k=0 

whereas the mobility operator can be written as 

^ oo 

M = Mo + Mo^-Zo- — -^gZoVfi, = f^o + Yl fj>'oVZo{-GZofgZoVti,, (6.32) 

where ^ 

Mo = ^ (6.33) 

is the one particle mobility matrix whereas Zo is one-particle convective friction matrix, 



given by the relation analogous to fl6.30p 

Zo = Zo- ZoV[iJ>Z,. (6.34) 

Since, similarly to the case of the Z operator, the velocity fields produced by GZo are 
force-free and torque-free, we obtain the relation 

ZoV = VZo = (6.35) 

which will be used in the following. 

Note that the scattering expansion fl6.32p is of the form (15. 6p with = HgZo, 
My = Zof-ig, Mo = ^Iq, and M = —Zo. Analogous scattering expansions for the 
kernels C and C introduced above read [25] 

oo 

C = ZoVfio - ZgZoViio = J2{-^oG)''2:oVfio, (6.36) 

oo 

C = iiJ>Zo - [loVZoQZ = ^ fi^VZoi-gZo)". (6.37) 

To obtain the response of the system on a macroscopic level, we need to average the 
above-defined hydrodynamic kernels over an ensemble of particle configurations. Next, 
the reduction procedure is carried out, according to the method outlined in Section [5l 
The kernels are reduced analogously to A in Eqs. (15. 7115. 9p . Using the scattering 
expansions ( ]6.10p .( !6732l) .( ]6.36p .( ]6.37p one obtains 

< m" >=< m" >'"■ + < G < C* > (6.38) 

< C >=< C >'''^ -<C G<Z> (6.39) 

< C >=< C - <Z G<C> (6.40) 
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and 

< Z >=< Z >^''^ -<Z G<Z> (6.41) 

These relations may be used to transform the response equations introduced in the 
previous section. For example, if the constant force E is applied to the particles, by 
averaging Eq. (16.281) one gets 



< / >=< > E- < Z > = 

(< (5* -<Z G<C' >)E - (< J -<Z >"■'■ G<Z >)vo 



(6.42) 



The above may be written in the form 

< / >=< C* E-<Z V > (6.43) 

where < v{r) > is the average velocity of the suspension as a whole 

<v>=vo + G<f>. (6.44) 

As it is seen from (16. 1116. 4p . the suspension velocity field v{r) has a simple interpretation: 
it is equal to the fluid velocity if r is inside the fluid and coincides with the rigid body 
motion wherever r lies inside the particle. 

Another quantity of interest is the average velocity of suspended particles 

U = l<5^C/.> (6.45) 

i 

which can be obtained by averaging Eq. (16.241) . using the reduction formulae (16.381) 
and (I6.39P and introducing the average suspension velocity according to (I6.44p . Such a 
procedure leads to: 

< U >= ^(< ^+ < E^* ^ >)' (6.46) 

7. Transport coefRcients 
7.1. Sedimentation and diffusion 

One of the fundamental problems in the physics of suspensions is the sedimentation 
phenomena - i.e. response of a suspension to a force field, e.g., gravity. The basic 
quantity here is the sedimentation velocity coefficient K, the ratio of the average particle 
velocity U to the acceleration of the external force field, E 

A- 4 (7.1) 

It is important to note that the sedimentation velocity is measured in the reference 
frame in which the fiuid as a whole is resting, i.e. < v >= 0. In this case Eq. (16.460 
gives 
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For the isotropic system, < ^ /x** >"''" is proportional to the unit tensor and the 
sedimentation coefficient may be then expressed as 

^=^Tr<5^Mg>'- (7.2) 

Moreover, this allows one also to find the collective diffusion coefficient, which is 
connected to K by the relation [28j 

O. ^ If*- (7.3, 

7.2. Viscosity 

The effective viscosity of a suspension, rjeff is obtained from the relation between the 
average stress of the system and the average rate of strain 

a = 2r^Xg,ff (7.4) 

with the effective value of the strain, ge//, given by 

9eff = ^['^a<V>/3+V,3<V>a] (7.5) 

In Eq. (17.41) . the tensor X is the fourth rank isotropic tensor, traceless and symmetric 
in its first and last index pairs: 

1 2 
2r = -(5a^5/3i. + - -Sap6f,^) (7.6) 

The stress in the suspension has two components - from the fiuid itself and from 
the force densities on particle surfaces [29], i.e. 

a = o-^^™'^ + o-P""* (7.7) 

with the particle contribution given by the ensemble average of the stresslet 

^part J2 SiS{r - Ri) > (7.8) 

i 

The partition (17.71) allows one to write the effective viscosity in the form 

Veff = ri + /\ri 

To calculate the effective viscosity, let us consider a problem of finding the force 
density / for the given fiow Vq in the absence of forces, .T-" = 0. This is a special case 
of (ESHD leading to 

/ = -Zv,. (7.9) 
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In particular, in the viscosity problem, one considers a linear velocity field of the form 

Vo = g-r (7.10) 

(with a symmetric and traceless matrix g) and looks for the stresslet, Sj of the induced 
force The response equation linking the local values of gj with the induced stresslet 

S. = 5^<g, (7.11) 

j 

defines the operator 

^ -pdzvj (7.12) 

where the projection operator defined in (I6.2ip has been used. 

The next step is to take the average over the particle configurations. Eq. (16.431) 
gives then 

<f>=-<Z>vo = -<Z >^"'^< V > (7.13) 

The stresslet may be obtained by acting on the above with the projection operator Pf. 
Expanding the fiow field in gradients and taking the lowest term leads to the following 
relation between stress and strain as 

a^-* = 1 < 5^ /xjf ge// (7.14) 

where the relation (I7.12p has been used. 

For the isotropic system the average tensor < yu'^'^ must be proportional to I, 

thus 



8. Fourier space formulation 
8.1. Sedimentation coefficient 

The transport coefficients defined above are often calculated using Fourier transform. 
In the case of the sedimentation coefficient, one starts with the Fourier transform of 
Eq. fl6.24p . which in the absence of an external fiow reads 

<V{k) >= K{k)E{k), (8.1) 

where 

^ik) = -J2^,^''''^ (8.2) 

i 

and 

K{k) = k- < ti'\k) > k (8.3) 
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is the wavevector-dependent sedimentation coefficient. In the above, 

^^'\k) = j^J2^'¥'^'^' (8-4) 

ij 

The usual sedimentation coefficient is then obtained as ^ hmit of (18. 3p 

K = -Tr hm < /x"(fc) > (8.5) 

3 ^0 

It is important to realize that the limit limfc^o < Ai**(fc) > in the above relation cannot 
be replaced by the A; = value of the kernel, < ^t**(A; = 0) >. This is caused by the 
presence of long-range hydrodynamic interactions in the system. Namely, the propagator 
G{r—r') contains terms which decay asymptotically as |r— r'|'>' with 7 < 3. While trying 
to calculate fc = value of the kernels, those long-range terms give rise to diverging 
integrals. 

An alternative way of calculating the sedimentation coefficient would be to start 
with the Fourier transform of Eq. (I6.46P 

< U(fc) >=< At"(fc) >''^'^ E+ < C\k) v{k) >, (8.6) 

with 

^(^) = ^E / C*(r)e^^-(^-^)dr (8.7) 



and then use the zero net flux condition [30] 

v{k = 0)=0 (8.8) 
which holds for incompressible fluid placed in an immobile container. This gives 

K = ^Tr lim < ^"(fc) >^"'^= ^Tr < ^"(fc = 0) (8.9) 

which is equivalent to (17.21) and does not involve small wavenumber limits, which makes 
it much more convenient in calculations. This time the value at = is well- defined 
since the long-range terms are absent in irreducible kernels [30] and thus those kernels 
are continuous at fc = 0. 

8.2. Viscosity 

The Fourier space formalism may be also used to define the viscosity coefficient. First, 
using the Fourier transform of the Oseen tensor 

G{k) = ^{l-kk), (8.10) 

one writes the velocity field in the absence of the particles as 

riPvoik) = il-kk)f,ik). (8.11) 
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The analogous relation between the average flow field in the presence of the particles, 
< i; >, and the external force density, /q will then define the wavevector dependent 
effective viscosity function rjeffik) 

Veff{k)e < v{k) >= (1 - fcfc)/o(fc). (8.12) 

Again, the hydrodynamic viscosity coefficient is defined as the long wavelength limit of 

Veffik) 

Veff = limr/e//(fc) (8.13) 

The function rjeff{k) may be expressed in terms of the hydrodynamic kernels defined 
above. To this end we note that the fiow field in the presence of particles may equally 
well be expressed as 

7]k^ < v{k) >= (1 - kk){< f{k) > +/o(fc)). (8.14) 

Inserting the Fourier transform of Eq. ( 16. 28^ yields (for the homogeneous system in the 
absence of external forces) 

7]k'^ < v{k) >= (1 - kk){fo{k)- < Z{k) > vo{k)). (8.15) 

In the above, the Fourier transform of the kernel Z is defined as 

Z{k) = j e-'^ '' Z[r - r') e^'^'^'dr. (8.16) 

where we used the fact that for a homogeneous system Zir^ r') = Z{r — r'). 
Next we eliminate t^o(fc), using the identity 

Vq = — ^ — <v> (8.17) 

1-G<Z> 

obtained by combining Eq. fl6.44p with Eq. (17. 9p . 
Finally 

+ (1 - kk) ^ < Z{k) >) < v{k) >= (1 - kk)f,{k). (8.18) 
1 - G < Z{k) > J 

Comparing Eq. fl8.12p with Eq. (I8.18P we obtain 

A77 = 77eff - 77 = lim 4tt ( (1 " ^fc) < Z{k) > ) : (1 - kk) (8.19) 

The above relation again involves a cumbersome fc ^ limit which cannot be replaced 
by the corresponding value at = 0, not only because of the l/k"^ term in fl8.19p but 
also since < ^ > is a long-range kernel, ill-defined at = 0. However, Eq. (I6.4ip gives 
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thus the relation fl8.19p may be rewritten in terms of the irreducible kernel < Z{k) >*^'' 

^■q = lim ^ ((1 - kk) < Z{k) : (1 - kk) (8.21) 

An explicit expression for the above limit may be obtained by expanding Eq. (18.161) 
in k and using the fact that the fields produced by the operator 2 are force free and 
torque-free. Thus the lowest order term in this expansion is 0{k'^) and corresponds 
to the third multipole (stress-strain) of force and velocity fields as defined in (16.131) 
and ( I6.14p . The coefficient in this term is thus proportional to the right hand side of 
Eq. (I7.12p . and the proportionality constant may be obtained by isotropy considerations 
(a detailed derivation may be found in Ref. [31]). Finally: 

= 1^ &o < 5: ^f^" "^' >^V= < E < >^^^V (8-22) 

where the last equality follows from the fact that irreducible kernels have a well-defined 
value at fc = 0. 

9. Linear response for Smoluchowski dynamics 

The above developed formalism may also be applied to the calculation of the system 
response in the long-time regime, when the memory effects (caused by the relaxation of 
distribution of particle positions) become important. As a first step towards solving this 
problem, we apply the linear response theory to generalized Smoluchowski equation, 
which governs the evolution of the particle distribution function in the configuration 
space, P{X,t) for the colloidal suspension. In the absence of external disturbances, the 
equilibrium distribution is given by 

Pe,(X) = e-'^*W/Q, (9.1) 

where is the potential of interparticle forces. 

Next, we disturb the system by introducing the imposed fiow field VQ{r) and 
external forces £ = {Ei,...,En) and calculate an induced mean force density 
and particle current. The evolution of P{X,t) is then given by the Generalized 
Smoluchowski Equation [28] 

^p(x,t) = r»(x,t)P(x,t) 

where the Smoluchowski operator, V{X,t), in the presence of the fiow Vo{r) and 
external forces S reads 



N o 



OR, 



+ ^-C/(X)-^o. (9.2) 
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Here D{X) is the diffusion matrix 

Di, = kBTfif^, (9.3) 

and 

Ff* = -V.<^ (9.4) 

are the interparticle forces. 

For later use, we also introduce the adjoint Smoluchowski operator, C, which obeys 

T>P,g{X)... = Peg{X)C... (9.5) 

Next, we find the mean particle current and force density. The former is given by 
the following ensemble average 

N N 

< j{r, X) >t=< R^Kr - R^) >=< Yl ^^^^("^ - >f (9-6) 

i=l i=l 

where the symbol < >t denotes the average over P{X,t). Inserting the explicit form of 
adjoint Smoluchowski operator yields 



< j{r, X) >,=< J2 (^''gx + + ■ ^^^^ + ^^"^ ~ (^-^^ 

i=i ^ 



where :F^"* = (F^;"*, F^"*, F*^*) and { }i stands for i-th component (in particle 
indexes) of the operator in brackets. For example 

{£ ■ /x(X)}, = J2Ej- Mi. = E ■ (9.8) 

j 3 

where the symmetry of mobility matrix has been used in the last equality. Moreover, in 
order to keep the notation simple, from now on we denote translational part of mobility 
matrix /x" simply by ^t, as only /j," appears in subsequent considerations. Analogous 
convention applies to C and C*, which will be written as C and C respectively. 

By considerations similar to the above one can also find the mean force density. As 
it has been shown in [32] it is given by the formula 



< /(r, X) >t=< (r'^ + J' + £)- C{X) - Z{X)vo >t . (9.9) 

In deriving the linear response formulas for the system of Brownian particles the 
approach due to Felderhof and Jones [331 El] is adopted. It is assumed that particles 
were at equilibrium in the infinite past so that 

P{X,t^-(x>)=P,,{X) 

Subsequently the fields £ and Vq are turned on and the distribution changes to 

P(X, t) = P,,{X) + 6P{X, t), (9.10) 
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with 6P{X,t) obeying (to the hnear order in £ and Vq): 

^^^^^ - T^^P = ■ [(m^W + Cvo{t))P.,] . (9.11) 

The solution of the above equation with initial condition 6P = for t = — oo is given 
by 

6PiX,t) = -P,,j' rft'e^(*-*')(_^ + ;3:F) . [^JiS{t')+Cv,{t')]. (9.12) 
This allows us to rewrite the expressions for < /(r, X) >t and < j(r, X) >t as 

< j{r) >t=< j >r* + <j >r*= j dr'{Y,Eir,r')E{r',t) + Y,,{r,r')vo{r' ,t)) + 

+ fdr' [ dt'{X,E{r, r', t - t')E{r', t') + X,„(r-, r', t - t')vo{r', t')), 

(9.13) 

< fir) >t=< f >r* + < / >r = / dr'{YfE{r,r')Eir',t)+Yf,ir,r')vo{r',t)) + 
f dr' I dt'{XfE{r,r',t~t')E{r',t') + Xf,{r,r',t-t')vo{r',t')), 

(9.14) 

where an auxiliary force field E{r,t) was introduced, such that 

E,{t)= [ 5ir-R.^E{r,t)dr. (9.15) 



and we have singled out instantaneous and retarded part of system's response 
(corresponding to averaging over Peg and SP in Eq. (19.101) . respectively). The former 
contribution appears immediately after E or Vq is turned on and follows the change of 
the external perturbation, while the latter describes memory effects due to the change 
of the distribution function induced by external forces. 

Instantaneous response kernels introduced above are defined as follows 



N 



YjE{r, r') =< J2 - Ri)fiijS{r' - Rj) >, (9.16a) 

N 

Yj,{r, r') =< J2 - R^)C^{r') >, (9.16b) 

N 

YfE{ry) =< Y,Cir)j6ir' - R,) >, (9.16c) 

Yf,{r, r') =< -Z{r, r') >, (9.16d) 
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whereas time-dependent response kernels X are given by 

Af 

X^E{r, r', t) = -p-' < ^ 5(r - R,)[fx- V]^e^*[(V ■ f^%S{r' - R,) >, 

(9.17a) 

Af 

Xj,{r,r',t) = -p-^ < 5^5(r - i?i)[^- V]ie^*(V +f3:F) ■ C{r') >, (9.17b) 

i=l 

N 

XfE{r,r',t) = -13-^ < C(r)- V e^*^[(V ■ lAj^W - Rj) >, (9.17c) 

Xf,{r, r',t) = -p-' < C(r)- V e^*(V ■ C(r') >, (9.17d) 

where the symbols V and V denote the operator d/dX acting to the left and to the 
right respectively. 

10. The reduction of response kernels 



The instantaneous response kernels defined in fl9.16ati9.16dp are reduced according to 



the general formula ( 15. 7p . Using the scattering expansions ( ]6.10p .( !6732|) .( ]6.36p .( ]6.37p 
one obtains then 

< Yab >=< Ya. + < Ya. >'^^ G < YfB >, (10.1) 

where A = j,E and B = E,v. 

A somewhat harder task is to perform the reduction of the retarded response kernels 
X given by Eq. (19.171) . The general form of those kernels is 

X =< Ae^^B > 

with two operators A and B on both sides of the evolution operator e'^*. It is precisely the 
presence of the evolution operator in the kernels that makes the reduction complicated. 
The procedure is as follows. First, the adjoint Smoluchowski operator 

c = V +:f] ■ /X- V (10.2) 

is decomposed as 

AT 

C{1, 2,...,N)=J2 ^o(^) + SC{1, 2,...,N), (10.3) 

where Co{i) is the single particle operator 

Co{t) = DoVl (10.4) 
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with Vf denoting the Laplacian with respect to Ri and Dq - single particle diffusion 
coefficient. It is worth noting that Cq does not introduce any correlation between the 
particles. Thus the evolution operator can be written as a series 



S{t)+ / dTS{t-T)5LS{r)+ dr dT'S{t-T)6LS{T-r')6LS{T) + ..., (10.5) 
Jo Jo Jo 



with 

N 

S{l,2,...,N;t) = l[S{t;t), (10.6) 

i 

and 

S{i;t)^e^''^^'. (10.7) 

Next, the scattering expansions of the operators A, B and 5L are performed. Then, 
after inserting the expansions into < Ae^^B > one ends up with the representation of 
the retarded response kernel as a sum of terms of the following structure 



dl...ds / dr dr2... drnA'{t)S{t - n) 

Jo Jo Jo (lO.i 

dL\n)S{n - T2)6L\t2) . . . S{Tn)B\Tn)c{l, . . . , s) , 



where A', B and SL' stand for some elements of the scattering expansions of A, B and 
5jC respectively and s is the number of particles appearing in the given term. The time 
variables have been added to time-independent operators 6L', A' and B' just to indicate 
their positions relative to the evolution operators in the above integral. 

11. Diagrammatic expansion for time-dependent kernels 

Since the scattering expansion of time-dependent response kernels involves more 
operators than the instantaneous response kernels, we need to introduce new elements 
into the diagrams, namely: 

• the single-particle evolution operators S{i, r — r') are represented in the diagrams 
by horizontal sohd lines (e-bonds): 



r' 



• a dagger line 



represents the two-body interparticle forces ( J- - bond) 
single arrows (^,^) represent the operators V and V respectively 
double arrows (=^,«^=) represent (3~^ V and (3~^ V respectively 

For example the diagram 
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1 

stands for the kernel 



T 



T 



b 







(D 5) 



dld2 /i(12)/i(34)i(l,2,3,4) dr dr' 



[11.1) 



S{1, 2, 3, 4; i - t)5£i(2, 3, 4)5(1, 2, 3, 4; r - r')SC2{l, 2)5(1, 2, 3, 4; 0^(1, 2, 3, 4), 
with the corresponding blocks given by 



i(l, 2, 3, 4) = - 5(r - Hi)M<(l)G(12)M(2)G(21)M(l) 

G(12)Z„(2)G(23)Z„(3)G(34)J„(4)G(43)M>(3) Vs, 

SCi{2, 3, 4) V4 M<(4)G(43)M(3)G(32)M(2)G(23)M>(3) Va, 

^2(1,2) =F(21)M(1)G(12)M(2)G(21)M>(1) Vi, 

5(1,2,3,4) =- V4 M<(4)G(43)M(3)G(34)M(4)G(42) 
M(2)G(21)M>(l)5(r - R^). 

The exact form of M, and iVf^ depends on the specific kernel to be represented 
by the diagrams. For example, in the case of XjE, we put 



and 



M(i) - Zo{i) 



M>{i) = Zo{i)v{i)fXo{i) 



(11.2) 
(11.3) 

(11.4) 



As it is seen, the scattering sequences in time-dependent diagrams have a more 
complicated structure than those encountered before, not only due to the presence of 
several independent blocks, but also due to the appearance of divergence operators V 
and V- 
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12. Reduction of time-dependent diagrams 

The next task is to perform the reduction of the time-dependent diagrams along similar 
lines to the approach presented previously - i.e. by identification of connection lines. 
The definition of a connection line is analogous to that in instantaneous response terms: 
an operator G{ik, ik+i) is called the connection line of a term Rs{A, c) if the latter can 
be written as 



J Ai(ii, ^2, ■ ■ ■ , ik)ci{ii, ^2, ■ ■ ■ , h)G{ik, ife+i)A2(«fc+i, ■ ■ ■ , «s)c2(«fc+i, ■ ■ ■ , is)dld2...ds, 

(12.1) 

so that after the removal of G(ik,ik+i) the term i?s(A,c) becomes a product of two 
independent integrals. Integrals over time have been omitted in the above expression as 
they arc irrelevant to our definition. The nodal line and nodal blocks for terms Rs{A, c) 
arc also defined analogously to the instantaneous response case. Thus, for example a 
diagram of the form 



(D 6) 



t 



T 



T 







has a single connection line (the one joining the particles 2 and 3). This is also a 
nodal line of this diagram. 

However, because of the fact that retarded response terms consist of a number of 
individual operators A', 5L'{ti), 5L'{t2) . . . the nodal structure of -Rs(A, c) is usually very 
complicated and in general it is impossible to apply the concept of block distribution 
function here. To analyse the nodal structure of retarded response kernels an ordering 
of the graph nodes is introduced first. Namely, moving along the graph from the left to 
the right we index all the nodes with the subsequent natural numbers. For the diagram 
(D 6) one gets 
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- - 



t3 



(D 7) 



t 



T 



T 







The above defined ordering allows us to introduce the notion of a proper diagram. 
To define it, let us consider a diagram Rg with a scattering structure A(ii, i2, ■ ■ ■ , is) and 
a nodal line G{ik,ik+i) such that 



Mil, 12, 



Ai(ii,i2, 



,ik)G{ik,ik+i)M{ik^ 



.is)- 



[12.2) 



The diagram Rs will be called proper if all the operators in which the particles from 
{ii,i2, ■ ■ ■ ,ik} appear have smaller indexes than these in which {ik+i, ■ ■ ■ ,is} appear. 
Thus the diagram (D 2) is not proper whereas the one of the form 




(D 8) 



t 



r 



T 







is proper. Note that the definition of a proper term concerns only the scattering 
structure in a diagram, the correlation structure is irrelevant here. 



12.1. Nodal structure of time dependent kernels 

As the evolution diagrams consist of many different building blocks {A, B and 5C) 
no wonder that their nodal structure is much more complicated than that of the 
instantaneous response diagrams analyzed in Section [51 For example the diagram 
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5 fl^ — -n- 



2 *Gj-^ 




<3 



T 



D 

t' 



(D 9) 



has the nodal structure of the form 





4,5 





In graph theory the above structure is called a tree: a connected graph which do 
not contain any circuits (the lack of circuits stems directly from the definition of the 
nodal line). Unfortunately, the presence of many branches makes it impossible to apply 
in this case the methods developed in Section [51 In particular, the block distribution 
function cannot be defined on the nodal structure like that of the diagram (D 4), as it 
lacks the linear ordering. 
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Luckily, the nodal structure graph of the proper terms is simpler. Namely, in the 
proper diagrams by the definition left-right ordering of the vertices is compatible with 
the nodal structure. Hence the nodal graph of a proper diagram is a simple chain - a 
tree with two terminal vertices only. 

For example, the proper diagram of the form 




has a chain-like nodal structure graph of the form 



5 



3,4 



1,2 



Thus in proper diagrams, nodal lines divide the particles ii, . . . ,is into nodal blocks 
Ci, 6*2, . . . which can be ordered according to the place in the chain. This means that the 
nodal structure can again be written in the form Ci\C2\.--\Ck, where Ci, C2, ■■■■,Ck come 
one after another in the time integral (110.81) . For such a structure a block distribution 
function can again be defined by Eq. (14. 2p . 

As it was mentioned, these concepts cannot be applied in the case of improper 
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terms. However, it may be shown [35J that in the thermodynamic hmit the sum of all 
improper diagrams in the expansion of a given time-dependent kernel vanishes. (We give 
see the sketch of the proof in the Appendix.) Therefore in the subsequent analysis we can 
consider proper terms only. The fact that the time-dependent diagrams have a chain-like 
structure is an important result, since it allows us to use a concept of block distribution 
function and carry out the reduction procedure in the case of time-dependent response. 
This element was missed by the authors of Ref. |Mj who applied directly the block- 
distribution function analysis in their studies on linear response theory of viscosity, 
without showing first that the structure of the terms in respective scattering expansion 
is indeed chain-like. 

Because of the chain-like form of the diagrams, it is now relatively easy to sum the 
proper terms which share a similar nodal structure. For example, the proper diagrams 
of the kernel 

X =< Ae^^B > 
may be divided in the following groups 

(i) Diagrams with the articulation line in A-block 

(ii) Diagrams with the articulation line in 6L -block 

(iii) Diagrams with the articulation line in B-block 

(iv) Irreducible diagrams. 

Below, the reduction procedure is carried out for the diagrams of each type 
(i) Proper diagrams with the articulation line inside A-block are of the form 



A' 



A> 5C 

Gi— articulation line 



6C 




where the ovals stand for correlation functions and the divergence operators in each 
block are marked 

The kernel A may be now reduced analogously to (15.71) which gives 

Xi =< A< G < A>e^'B > (12.3) 
(ii) The proper diagrams with the articulation line inside 5jC-hlock are of the form 
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Thus, after the reduction, the diagrams of X2 sum up to 



Xo 



dr < Ae^(*-")5L< >"'^ G < 5L>e^^B > 



:i2.4) 



(iii) Nonvanishing diagrams with articulation hne inside B-block are of the form 



A 




B> 

G—art. line 



and they sum up to 

X3 =< Ae^'6B> G<B> > . (12.5) 
(iv) Finally, the irreducible diagrams give 

X4 =< Ae^'5B . (12.6) 
Eventually, summing up fll2.3tll2.6l) we get for the kernel X the following expression 
X(t) =< Ae^'6B >=< Ae^'6B + < A< G < A>e^'B > + 

t n 2 7) 

dr < Ae^(*-")5£< G < 6C>e^^B > + < Ae^'B< G < B> > 

The above algorithm may be now used to reduce the time-dependent kernels defined 
in Eq. fl9.17p . Namely, the analysis of the scattering structures of both retarded and 
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instantaneous kernels fl9.16aH9.17dl) leads to 

XAB{t) = XYsit) + Y-ZGXfB{t) + f drXY^t - r)GX^^(r) + X^ZmYf^. 

(12.8) 

where again A = j,E and B = E,v. 
13. Effective equations 

The reductions of instantaneous kernels Y and time-dependent kernels X carried out 
above may now be used to obtain the effective equations governing the dynamics of 
suspensions. Namely, using (110. ip in Eq. (I9.13P we get the following expression for the 
instantaneous part of the current 

< j >r= Y,eE + Y,,Vo = Y%E + Yj:v^ + YJ^G < f >'^' (13.1) 

where the definition of the instantaneous force density (19.140 was used. 

The retarded part of the current may be similarly obtained from the reduction 
formulae ^M)- 

< 3 >r = f dt'(X;-(t - t')E{t') + X^;:{t - t')vo{t')) + 

J — oo 

+ YJ:G I dt'(XfE(t - t')E{t') + Xf,{t - t>o(t')) + 

+ j At' dTX'^:{t-t' -T)G{XfE{r)E{t') + Xf,{r)v^{t')Y 
+ f <lt'XX{t-t')G{Y fEEit') + Y f^v.it')). 

J — oo 

The third term can be simplified by first changing the variables of integration to 
(t',t" = t' + r), then changing the order of integration, and finally using the fact that 
(c/. Eq.EH 

f dt'(XfE{t" - f)E{t') + Xf,{t" - t')v,{t')) =< f >[,f . (13.3) 
By this means Eq. (113. 2p can be rewritten as 

<J>r= r dt'{XJ^{t - t')E{t') + X-^{t ~t')vo{t')) + 

J — oo 

yj:g < f >r* + /* dt"xi^{t - t")G < f >t. . 

J — oo 



13.4) 



The equations for the instantaneous and retarded current are then added to yield 
the total current. The structure of the equations can be most clearly seen after the 
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Fourier transform in time: 

< jm >=< jm >^"^* + <m {yye+xy^{u))e{uj)+{yj:+x-:{oo))v{uj) 

(13.5) 

where we used the fact that the total suspension velocity may be written as 

<v>=VQ + G<f> (13.6) 

The Fourier transforms in time introduced above are defined as 

1 /"°° 

v{uj) = — <v>t e^'^Mt, (13.7) 
and analogously for j{uj)^ whereas the kernels X =< Ae^^B > are transformed as 



oo 



X{lu)= < Ae'^'B > e^^'dt. (13.8) 







In an analogous way one may derive the equation for the average force density 
(cf.EUD, getting 



< fiu) >= {YY + XY)E{u) + {Y^ + X^) < v{u) > (13.9) 

The above result can be inserted into the Stokes equation to yield, after the Fourier 
transform in space, 



Vk' < v{k,u) >= {YY + XYe)E{u) + {Yj: + Xj:) < v{k,u) > (13.10) 
13.1. Small k expansions of response kernels 

In the long wave limit the tensor YJ^ takes a particularly simple form. Namely, the 
scattering expansion (I6.37P gives 

/N oo 
< J2^W[J2'^'^'o'PZo{-gZo)%{r') dr' (13.11) 
1=1 «=o 

However, the integral over r' may be replaced by the action of projection operator 7-**. 
Then, using Eq. fl6.35p . which implies that ZoT'^ = 0, we get 



N N 



Y,,{k = Oy"- =<J2SiR^)^^oi^)'Pi^)Zo{i)V\i) J2S{R,) >l=nl (13.12) 



1=1 i=l 



where the expression fl6.33p for one particle mobility matrix has been used. The next 
nonvanishing term in the expansion of Yjy{k) in k is the second order one 

Y.M'''" = nl + k^y, , + (13.13) 
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with the tensor y-^ of the form 

y^, = y\M + yW^-kk), (13.14) 

where yj^ and y*^, are scalars representing longitudinal and transverse part of t/^^, 
respectively. Next, since Y'^J^. is adjoint to YJ^, we get 

YfE{k = Of-"^ = Yj^{k = Oy"' = nl. (13.15) 

On the other hand, again using the property fl6.35p we get a simple result 

Yf^{k = Oy"- = (13.16) 

The small k expansions of operators Y je and Y read 

YjEiky^' = nl + k''yjE + ..., 

Yf,{ky'' = k^yf^ + ... (13.17) 

Analogous expansions are carried out for the time-dependent kernels 

X7;(fe, oj) = ex,,{u) + . . . , (13.18) 
X7;(fe, u) = k^x,,{uj) + . . . , (13.19) 
XfE{k,u:) = k''xfE{uj) + (13.20) 

Using these expansions in Eqs. (113.51) and (113. lOp one arrives at the following relations 
for the diffusion current and force density for small but finite A;: 

j(fc, oj) - nv{k, uj) = {yjE + XjEi^))E{k, tu) + P(l - kk) [y]^ + x)„{uj))V{k, uj), 

(13.21a) 

+ y). + ^).{^)>{k, uj) = {\- kk) (^/o(fc, u) + nE(k, u) 

+ k\y)E + x)E{'^))E{k, a;)) (13.21b) 
Here, yjE and XjE are defined by 

Y,E{k = = y,El. (13.22) 

and 

Xf^{k = 0,uj) = x,£;(cj)l (13.23) 

Moreover, y^^^ and x* denote the transverse part of the operators y^^i, and Xab respectively 
and the incompressibility condition (1 — kk)v{k) = v{k) was used. 

The dynamics described by Eqs. fll3.21ap and fll3.21bp is relatively complex. First 
of all, there are direct effects. First, an external force E applied to the particles induces 
the diffusion current 

j'^ = j-nv (13.24) 
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which is the particle current measured relative to the average suspension velocity frame. 
The intensity of that effect is measured by the sedimentation coefficient, which now 
becomes frequency-dependent and reads 

Kiu;) = ^iyjE + x,E{to)) (13.25) 

Moreover, as seen in (113.21bl) . the suspension velocity field is induced by the overall 
external force acting on the particles and the fluid 

Ft,t = fo + nE (13.26) 

The effective viscosity of the suspension is modified by the presence of the particles and 
reads 

r^^ff{uj) = 7^ + y}, + x%iuj). (13.27) 

Finally, there are also cross effects linking the suspension velocity v{uj) with the external 
force acting on the particles E{uj) and the diffusion current j{uj) with v{uj). The 
intensity of those couplings is measured by the coefficients y^j-^ + x^j^{uj) and ?/*^ + x*j,(ti;) 
respectively. However since f'E kernels are adjoint to jv ones (cf. Eqs. fl9.16dl) and 
(I9.17P ) the above coefficients are in fact equal. This is a manifestation of the Onsager 
symmetry as suggested by Nozieres [5^ . 

14. Summary 

The response of a composite system with field-induced forces was studied using a newly 
developed diagrammatic method. The method may be used in both instantaneous and 
retarded response analysis. It was shown that in both cases it is possible to describe 
the system's response by a set of transport coefficients which depend solely on local 
properties of the medium. The expressions for the transport coefficients obtained with 
use of the diagrammatic technique were shown to be well-behaved and free of divergences 
even in the presence of long-ranged forces. Thus they represent a proper starting point 
for calculation of the transport coefficients and for construction of approximate methods. 

As mentioned in the Introduction, a subsequent article [13] will discuss the 
application of the above methodology to the problem of the settling velocity and its 
fluctuations in a non-Brownian suspension. This task is more complex than analogous 
analysis for the Brownian suspension, presented in Sec. I7.H since the distribution 
functions in that case correspond to the nonequilibrium (though stationary) state. Using 
the diagrammatic technique one can derive correlation functions in this state. Again, 
the crucial element of the derivation is the reduction procedure. Due to its complexity, 
this procedure is nearly impossible to carry out were it not for the rigorous methodology 
provided by the diagrammatic method. 
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Appendix A. Simplification of time-dependent diagrams 

In this appendix we sketch the idea of the proof that the sum of all time-dependent 
improper diagrams vanishes. The detailed proofs may be found in [35] . 

First, let us introduce a few additional definitions concerning structure of the 
diagrams from the expansion of < Ae^^B >. First, let us note that when one removes 
all correlation functions and e-bonds from a given diagram, it decomposes into a number 
of subdiagrams - scattering blocks, representing A, B or subsequent 5C^s operators. 
The vertex in a given block which is most to the left (right) will be called first (or 
last) vertex of the block respectively. Finally, right (left) terminal block is a block with 
the property that the particle fine passing through its last(first) vertex v does not pass 
through any other vertex in a diagram more to the right (left) than v. 

For example, the diagram (D 6) consists of four scattering blocks. The first one 
(from the left) is a left terminal A block, the next one is a left terminal 5C block. Then 
there is another 5C block and finally a right terminal B block. 

Note that every improper diagram must contain one of the following: either a 
right (or left) terminal 5C block or a right terminal A block or a left terminal B block. 
Next, we consider these cases in order. 

The case when a diagram contains a right terminal 5C block is relatively 
straightforward. It suffices to note that every 5C block ends with the operator. 
In the case of the right terminal 5C block this divergence operator has nothing to act 
on to its right and thus the value of such a diagram vanishes. 

The case of left terminal 5C block is a bit more complicated. Let us denote such a 
block by Li,. There are two possibilities: 

a) Lb begins with o operator i.e. Vi A*o(^)'^(0-^o(0 

b) Lb begins with Fjifi^{i)'P{i)Zo{i) 

Here i denotes the particle with which Lf, begins whereas j is some particle from the 
diagram different from i. For the diagrams in (a) , using integration by parts one can 
transform o operator at the beginning of Lf, for — 1 ■ o operator. But, as Lf, is 
the left terminal block, after such operation, the differentiation in — l-O acts only on 
the correlation function on the far left of the diagram. However, for the equilibrium 
distribution functions 

-r'VX'll, . . . , s) = - ^ F,, <"(!, ...,s)- f F,(,+i)<''(l, . . . , s + l)d{s + l) (A.l) 

i=i 

Thus each diagram in (a) may be written as a sum of a number of diagrams in (b) with 
the same scattering structure, taken with an opposite sign. In this way one can show 
that the total sum of all diagrams in (a) and (b) vanishes. 

It remains to consider two more cases: the diagrams with a right terminal A-block 
and those with left terminal i?-block. However, in the first case, it suffices to transform 
O operator at the end of A-block to — 1 ■ O using integration by parts, and then note that 
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again the divergence operator has nothing to act on to its right. Eventually, in case of 
diagrams with left terminal B-block, the proof is analogous to the one concerning left 
terminal SC blocks presented above 

References 

[1] S. Torquato. Random Heterogeneous Materials: Micro structure and Macroscopic Properties,. 

Springer, Berlin, 2002. 
[2] M. Sahimi. Heterogeneous Materials. Springer, New York, 2003. 

[3] R. Landauer. Electrical conductivity in inhomogeneous media. In J.C. Garland and D.B. Tanner, 

editors. Electrical Transport and Optical Properties of Inhomogeneous Media, AIP Conf. Proc. 

No. 40, pages 2-43. AIP New York, 1978. 
[4] J. G. Kirkwood. On the theory of dielectric polarization. J. Chem. Phys., 4:592-601, 1936. 
[5] J. D. Ramshaw. Existence of the dielectric constant in nonpolar fluids. Physica, 62:1-16, 1972. 
[6] D. Bedeaux and P. Mazur. On the critical behaviour of the dielectric constant for a nonpolar 

fluid. Physica, 67:23-54, 1973. 
[7] B. U. Felderhof, G. W. Ford, and E. G. D. Cohen. Cluster expansion for the dielectric constant 

of a polarizablc suspension. J. Stat. Phys., 28:135 164, 1982. 
[8] C. W. J. Beenaker. The effective viscosity of a concentrated suspension of spheres (and its relation 

to diffusion). Physica A, 128:48-81, 1984. 
[9] B. Cichocki and B. U. Felderhof. Renormalized cluster expansion for multiple scattering in 

disordered systems. J. Stat. Phys., 51:57-76, 1988. 
[10] R. Balescu. Equilibrium and nonequilibrium statistical mechanics. John Wiley and Sons, New 

York, 1975. 

[11] R. G. Barrera, G. Monsivais, W. L. Mochan, and E. Anda. Diagrammatic approach to the effective 
dielectric response of composites. Physical Review B, 39:9998-10008, May 1989. 

[12] R. G. Barrera, C. Noguez, and E. V. Anda. A new diagrammatic summation for the effective 
dielectric response of composites. J. Chem. Phys., 96:1574-1581, January 1992. 

[13] B. Cichocki and K. Sadlej. (to be published). 

[14] R. E. Caflisch and J. H. C. Luke. Variance in the sedimentation speed of a suspension. Phys. 
Fluids, 28:259, 1985. 

[15] H. Nicolai and E. Guazzelli. Effect of the vessel size on the hydrodynamic diffusion of sedimenting 

spheres. Physics of Fluids, 7(1) :3 5, 1995. 
[16] H. Nicolai, B. Herzhaft, E. J. Hinch, L. Oger, and E. Guazzelli. Particle velocity fluctuations and 

hydrodynamic self-diffusion of sedimenting non-brownian spheres. Physics of Fluids, 7(l):12-23, 

1995. 

[17] P. N. Segre, E. Herbolzheimer, and P. M. Chaikin. Long-range correlations in sedimentation. 

Physical Review Letters, 79(13):2574-2577, 1997. 
[18] G. K. Batchelor. Sedimentation in a dilute polydisperse system of interacting spheres. Part 1. 

General theory. J. Fluid Mech., 119:379 408, 1982. 
[19] G. K. Batchelor and C. S. Wen. Sedimentation in a dilute polydisperse system of interacting 

spheres. Part 2. Numerical results. J. Fluid Mech., 124:495-528, 1982. 
[20] F. Feuillebois. Multiphase Science and Technology, page 763, 1991. 

[21] J. P. Hansen and I. R. McDonald. Theory of Simple Liquids. Academic Press, 1976. Sec. 4.6. 

[22] M. A. J. Michels. The convergence of integral expressions for the effective properties of 
heterogeneous media. Physica A, 157:377-381, 1989. 

[23] G. E. Uhlenbcck and G. W. Ford. The theory of linear graphs with applications to the theory 
of the virial development of the properties of gases. In J de Boer and G E Uhlenbeck, editors. 
Studies in Statistical Mechanics, volume 1, pages 119-211. North-Holland, Amsterdam, 1962. 



Diagrammatic approach to response problems in composite systems 42 

[24] P. Mazur and D. Bedeaux. Generalization of Faxon's theorem to nonsteady motion of a sphere 

through an incompressible fluid in arbitrary flow. Physica, 76:235-246, 1974. 
[25] B. U. Felderhof. Many-body hydrodynamic interactions in suspensions. P%szca, 151 A: 1-16, 1988. 
[26] B. Cichocki, B. U. Felderhof, and R. Schmitz. Hydrodynamic interactions between two spherical 

particles. Physico Chemical Hydrodynamics, 10:383-403, 1988. 
[27] R. Schmitz and B. U. Felderhof. Mobility matrix for two spherical particles with hydrodynamic 

interaction. Physica A, 116:163-177, 1982. 
[28] P. N. Pusey. Colloidal suspensions. In J. P. Hansen, D. Levesquc, and .J. Zinn- Justin, editors, 

Liquids, Freezing and Glass Transition, pages 763-942. Elsevier. Amsterdam, 1991. 
[29] S. Kim and S.J. Karilla. Microhydrodynamics. Butterworth-Heinemann, Boston, 1991. 
[30] P. Szymczak and B. Cichocki. Memory function for collective diffusion of interacting Brownian 

particles. J. Chem. Phys., 121:3329-3346, 2004. 
[31] B. Cichocki, B. U. Felderhof, and R. Schmitz. The effective viscosity of suspensions and emulsions 

of spherical particles. Physica, 154A:233-256, 1989. 
[32] B. U. Felderhof. Brownian motion and creeping flow on the Smoluchowski time scale. Physica A, 

147:203-218, 1987. 

[33] B. U. Felderhof and R. B. Jones. Linear response theory of sedimentation and diffusion in a 

suspension of spherical particles. Physica A, 119:591-608, 1983. 
[34] B. U. Felderhof and R. B. Jones. Linear response theory of the viscosity of suspensions of spherical 

Brownian particles. Physica, 146A:417-432, 1987. 
[35] P. Szymczak. Memory function for collective diffusion of interacting Brownian particles. PhD 

thesis, Warsaw University, Institute of Theoretical Physics, June 2001. 
[36] P. Nozicres. A local coupling between sedimentation and convection: application to the Beenaker- 

Mazur effect. Physica A, 147:219-237, 1987. 



